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ABSTRACT 

The orbital evolution and stability of planetary systems with interaction 
from the belts is studied using the standard phase-plane analysis. In addition to 
the fixed point which corresponds to the Keplerian orbit, there are other fixed 
points around the inner and outer edges of the belt. Our results show that for 
the planets, the probability to move stably around the inner edge is larger than 
the one to move around the outer edge. It is also interesting that there is a limit 
cycle of semi-attractor for a particular case. Applying our results to the Solar 
System, we find that our results could provide a natural mechanism to do the 
orbit rearrangement for the larger Kuiper Belt Objects and thus successfully 
explain the absence of these objects beyond 50 AU. 

Keywords: Phase-plane analysis; bifurcation; planetary systems; stellar dynamics. 
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1. Introduction 

Newton's Law of gravity and motion has been used as the most fundamental law for 
the Physical Sciences since its success in explaining the motion of celestial bodies in the 
Solar System. Thus, the Newton's Law was in fact first proved in the astronomical context. 
It was then apphed to other fields successfully. 

For example, Clausen et al. [1998] studied periodic modes of motion of a few body 
system of magnetic holes, both experimentally and numerically. Kaulakys et al. [1999] 
showed that a system of many bodies moving with friction can experience a transition to 
chaotic behavior. Xia [1993] did some work on Arnold diffusion in the elliptic restricted 
three-body problems. 

Many new discoveries of extrasolar planets have been made recently (Marcy ct al. 
[1997]) and these events really provide exciting and important opportunities to understand 
the formation and evolution of planetary systems, including the Solar System. According 
to the Extrasolar Planets Catalog (http://cfa-www.harvard.edu/planets/catalog.html) 
maintained by Jean Schneider, 117 planets have been detected till date (July 2003). These 
planets have masses ranging from 0.16 to 17 Jupiter mass (Mj) and semimajor axes ranging 
from 0.04 AU to 4.5 AU and a wide range of eccentricities. Therefore, the dynamical 
properties of some of these planets are very different from the planets in the Solar System. 

Nevertheless, there do exist similarities between the extrasolar and the Solar System 
planets. For example, there is a new discovery of a Jupiter-like orbit, i.e. a Jupiter-mass 
planet on a circular long-period orbit (Marcy et al. [2002]). Moreover, some planetary 
systems are claimed to have discs of dust and they are regarded to be young analogues of 
the Kuiper Belt in our Solar System (Jayawardhana et al. [2000]). 

If these discs are massive enough, they should play important roles in the origin of 
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planets' orbital elements. In fact, Terquem & Papaloizou [2002] provide an important 
mechanism to explain the observed planetary systems: massive planets form at the inner 
part of the system and move on a circular orbit due to tidal interaction with the disc, 
whereas lower mass planets form at the outer part and interact with the inner ones. 
Thus, the eccentricities of the massive inner planets increase, which is in accord with 
the observations. Moreover, the simulations in Jiang & Ip [2001] show that the orbital 
configuration of the planetary system of upsilon Andromedae might initially be caused 
by the disc interaction. Although the discs might be depleted and gradually become less 
massive, there could be belts surviving until now as the Asteroid and Kuiper Belts in the 
solar system and these belts would be important for whole dynamical histories of planets. 

On the other hand, Thommes, Duncan & Levison [1999] claimed that Neptune was 
closer to Jupiter and got scattered outward. Yeh & Jiang [2001] studied the orbital 
migration of scattered planets and completely classified the parameter space and the 
solutions. They concluded that the eccentricity always increases if the planet which moves 
on a circular orbit initially, is scattered to migrate outward. Thus, the orbital circularization 
must be important for the scattered planets if they are now moving on nearly circular 
orbits. 

Since belts of planetesimals often exist within a planetary system and provide the 
possible mechanism of orbital circularization, it is important to understand the solutions of 
dynamical systems with planet-belt interaction. Therefore, Jiang & Yeh [2003] did some 
analysis of the solutions for dynamical systems of planet-belt interaction for a belt with 
simple density profile c/r, where c is a constant. However, because the gravitational force 
from the belt involves complicated elliptic integrals and is difficult for the construction 
of theorems, they used simplified analytic formulas to estimate this force. Their main 
conclusion is that the inner edge of the belt is more stable than the outer edge of the belt 
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for the planets. Their result is actually consistent with the observational picture of the 
Asteroid Belt. 

In this paper, we complete the work initiated in our earlier paper Jiang & Yeh [2003], 
by numerically evaluating elliptic integrals to obtain the gravitational force of the belt 
without introducing any simplifications. Thus, this study is intended to improve the 
understanding about orbital evolution of a given planet-belt system and hopefully explain 
why there are less Kuiper Belt Objects larger than 160 km in diameter beyond 50 AU in 
the outer Solar System. In Section 2, we mention our basic governing equations and in 
Section 3, we study the locations of the fixed points. The phase-plane analysis is discussed 
in Section 4 and Section 5 concludes the paper. 



2. The Model 

The belt is an annulus with inner radius ri and outer radius r2, where ri and r2 are 
assumed to be constants. According to Lizano & Shu [1989], the belt's surface density scales 
with distance as r~^. We thus assume the density profile of the belt to be p(r) — c/r^, 
where c is a constant, completely determined by the total mass of the belt. One might 
notice that there are discontinuities at r = ri and r = r2. This is probably acceptable 
because the edges of the belt are rather sharp as one can see from a plot of the Asteroid 
Belt in Murray & Dermott [1999]. 

The units of length and time are chosen so as to easily compare with our Solar System. 
There are two belts for our Solar System: the Asteroid Belt in the inner part of the Solar 
System and the Kuiper Belt in the outer part. We set — 3 and r2 — 6 because (1) when 
the length unit is 1 AU (and the time unit is l/(27r) years), the interval [ri,r2] covers the 
region of the Asteroid Belt approximately, (2) when the length unit is 10 AU (and the time 
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unit is 10^/^/(27r) years), the interval [ri,r2] covers the region of the Kuiper Belt. 

We assume the distance between the central star and the planet to be r, where r is a 
function of time. The total force / acting on the planet, includes the contributions from the 
central star and the belt. The gravitational force from the central star is 

/. = (1) 

where we have set the mass of the central star to be 1 and also the gravitational constant 
G = l. 

On the other hand, the gravitational force from the belt for the planet is complicated 
and involves elliptic integrals. We use the equations in Appendix A to get the force 
numerically. 

Because there might be some scattering between the planetesimals in the belt and 
the planet, the planet should experience the frictional force when it goes into the belt 
region. This frictional force should be proportional to the surface density of the belt at 
the location of the planet, and the velocity of the planet. However, if the planet is doing 
circular motion, the probability of close encounter between the planet and the planetesimals 
in the belt is very small and can be neglected here. Thus, we can assume that the frictional 
force is proportional to the radial velocity of the planet dr/dt only, and ignore the d9/dt 
dependence. 

Therefore, the frictional force fa is proportional to the surface density of the belt and 
radial velocity dr/dt. Hence, we write down the formula for frictional force as 

~ dv 

fa = -ap{r)-, (2) 

where a is a frictional parameter and p is the density profile of belt. 

Because the dO/dt component of the planetary velocity is ignored in the frictional force. 
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the angular momentum I is conserved here, so we have 

mr^de = Idt. (3) 

This imphes that Kepler's second law is still valid here and, because of this, we can use 9 
as our independent variable. We use 9 to label time t from now on and one can easily get t 
from the above equation. 

Since u — l/r and using Equation (3), we have 

dr dr dd I dr I du 

dt ^ d9dt ^ mr^d9 ^ ~md9' ^ ^ 

Equation (2) can be rewritten as 



fa = -ap{r) 



I du 
m d9 



alp du 
m d9 



(5) 



Since the planet feels the frictional force only when it enters the belt region, we define 



fa Pfai 



(6) 



where 



1 if ri < r < r2, 

if r > r2 or r < ri. 



In general, the equation of motion for the planet is (Goldstein [1980]) 



fu 
d9'' 



m ] 



(7) 



where u = 1/r and m, I are the mass and the angular momentum of the planet and / is the 
total force acting on the planet. This equation is only valid when there is no non-radial 
force. We use the polar coordinates (r, 9) to describe the location of the planet. 
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Because f — fs + fb + fa, transform Equation (7) into 

(fu rn^fb (3apdu 

Therefore, the equation of motion for this system can be written as 



-u- Pki-^v + k2-k2^ = g2{u,v), (9) 

dO 



where ki = a/l, k2 = rri^/f 



3. Fixed Points 

The fixed points {u, v) of Equation (9) satisfy the following equations : 

-u - PkiA,v + k2- k2^ = 0. (10) 
From Equation (10), the fixed points {u, 0) satisfy 

- k2u'' + k2fb = 0. (11) 

The fixed points can be numerically determined for any given k2 and fb and it is 
obvious that the location of fixed points does not depend on the friction. Since at fixed 
points V — du/d6 — 0, the fixed points correspond to circular orbits. Particularly, when 
/b = 0, by Equation (11), the fixed point satisfies u = k2 and this actually corresponds to 
the Keplerian orbit. (One can check easily by writing /c2 as a function of usual angular 
velocity.) 
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The density profile of the belt is p{r) — c/r'^, where c is determined by the total mass 
of the belt (Lizano & Shu [1989]). Hence, the total mass of the belt is 

Mb = r p{r')r'dr'd(f) = 27rc(lnr2 - Inn). (12) 

Jo Jri 

Prom Appendix A, we can then calculate for any given total mass of the belt Mf, 
as the constant c can be uniquely determined from Mh. Moreover, since m is always 0.001 
(about Mj) in this paper, k2 is directly related to the angular momentum of the planet. 
Therefore, we use and k2 as the two main parameters to determine the number and 
locations of fixed points in this system. 

Figure 1 is the bifurcation diagram of fixed points and these curves show the locations 
of fixed points m as a function of k2 for a given M^. From this plot, we know that when the 
total mass of the belt is very small, = 0.01, the curve on the k2 — u plane is almost a 
straight hne with slope 1. That is, there is one and only one fixed point u k2. This can be 
easily verified by Equation (11) with small /b, and it corresponds to the Keplerian orbits. 
When the mass of the belt is not negligible, the curves form two needle-like structures 
pointing along both, the inner edge of the belt, where u— 1/ri = 1/3, and the outer edge of 
the belt, where ii = l/r2 = l/6, and there would be more than one fixed point for particular 
values of ^2- In addition to the one corresponding to the Keplerian orbit, other fixed points 
are located around inner and outer edges of the belt. Thus, the belt could bring the planets 
to non-Keplerian circular orbits by the gravitational and frictional forces. 

When Mft = 0.1, the needle-like structures imply the existence of two fixed points for 
0.5 < k2 < 0.7. When M}, — 0.5, the two needle-like structures become even bigger and 
there are three fixed points for 0.5 < A;2 < 1 , two fixed points around k2 — 0.06 and one 
fixed point for most other values of k2. When M5 = 1, the structure grows even larger and 
becomes completely different from a straight line. There are three fixed points for larger /c2 
(/c2 > 0.6) and one fixed point only for both intermediate k2 (0.1 < k2 < 0.5) and tiny k2 
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(/c2 < 0.02). Since the needle- like structure around the outer edge is tiny, the probability 
that the planet moves stably around the outer edge is much smaller than near the inner 
edge. 



4. Phase-Plane Analysis 

Following the linearization analysis, the eigenvalues A corresponding to the fixed points 
(it*, 0) satisfy the following equation: 



'dgi 
du 



(«.,0) ) \ov 



-A - 



(«*,o) 







(dg2 




\ dv 


K,o)) 


\ du 





0, 



(13) 



where gi{u,v) — v and g2{u,v) = —u — f3{kipv)/u^ + k2 — {k2fb)/u^- Hence, the above 
equation becomes 



A^ + ^A+(l-^ 



, dA 

{u„0) OU 



(W.,0)y 



= 0, 



where A — 



When the planet is out of the belt region, that is P — 0, we have 

A = ±. 

and there are two possible cases: 



'1 



-1 - fca^ 

OU 



(«*,0) 



Case I: If -1 - fcaf^ 



Case II: If -1 - fcafj 



(w*,0) 



> 0, this fixed point is a saddle point. 



(«*,o) 



< 0, this fixed point is a center point. 



When the planet stays in the belt, that is /9 = 1, we have 



2ul I \ 



Ui 



(14) 
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We define A = 1 - 4 1 + A;2^^ 



du 



V^j ' ^^^^ ^® have three possible cases: 



(u*,0), 



Case 1: When A < 0, two eigenvalues Ai^2 — —a±bi for a, 6 > 0, so the fixed point is a 
stable spiral point. 

Case 2: When < A < 1, both eigenvalues are negative, so the fixed point is a stable 
node. 

Case 3: When A > 1, one eigenvalue is negative and the other is positive, so the fixed 
point is a saddle point. 

We understand that when the real parts of the eigenvalues of a fixed point equal to 
zero in the first-order linearization analysis, the properties of this fixed point should be 
determined by the higher order analysis in general. However, to simplify our language, we 
still use the term "center point" for any fixed point having an eigenvalue with zero real 
part. This is a good choice because our numerical results indicate that these points are in 
fact the center points. 

We now derive the mathematical formula for A. At first, from Equation (14), we 
calculate (please see details in Appendix B) 



OA 
du 



{u„0) 



ul du 



_2h 

(«.,o) ul 



/ r'p(r') 



uj dr 
1 + r''^ul 



-E 



1 



—F\dr\ 



(15) 



where E is an elliptic integral of the second kind and F is an elliptic integral of the first 
kind. Hence, we have 



ul V . dA 



(w*,0)y 



\kxp) I Jr> 



1 + /•'-(/ 



./2 ,,2 



-E 



1 



= 1 



+ 



^u^ 



r p{r' 



(1 -r''u*)2(l + r'ii*) l + r'ii* 

1 + ,-,,= ^ 1 



dr. (16) 
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If we calculate the value of A numerically, we can then understand the properties 
of this fixed point. For example, in Figure 2-6, we choose four sets of (Afj, A;2) values 
and make the phase-plane plots, orbits etc. to demonstrate the orbital evolution, and the 
stability of the system to provide a general picture for the readers. In the results below, 
the frictional parameter a were chosen to be small numbers but large enough to make the 
solution curve not too tight to be studied. The value of the frictional parameter a is related 
to the viscosity of the belt and is too complicated to be well-determined observationally. 
The values of a in the models would affect the time scales of orbital evolution but would 
not change the general understanding of the dynamical properties of the system. For the 
parameters we choose, some results show that the orbits become nearly circular in an order 
of 10^ years if the unit of length is 10 AU. This is rather short and perhaps only possible 
during the period in which the planetary system is just formed. Nevertheless, one can easily 
choose smaller values of a to get results of orbital evolution with longer time scales when 
the results are used to describe more recent systems. 

Figure 2 is the case when we set M5 = 0.5 and k2 = 0.6. We assume the frictional 
parameter a = 0.001 for this case. We know from Figure 1, that there must be three fixed 
points and Figure 2(a) shows that the one around u — 0.6 is a center point and the other 
two near the inner edge of the belt are spiral and saddle points, respectively. The blue 
dotted curve around the center point in Figure 2(a) is then plotted on the x — y plane as the 
blue dotted orbit in Figure 2(b). It is a precessing elliptical orbit with small eccentricity. 
The red solid curve approaching the spiral point in Figure 2(a) is then plotted on the x — y 
plane as the red solid orbit in Figure 2(b). This orbit keeps changing the semimajor axis 
and eccentricity and finally becomes a circular orbit with radius r — ri, i.e., close to the 
inner edge of the belt. Figures 2(c)-(d) show u and v as functions of 9 for both above orbits. 
Because u becomes a constant and v approaches to zero when ^ > 30 for the red solid 
curves in Figures 2(c)-(d), these two figures reconfirm that the red solid curve in Figure 
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2(b) finally becomes a circular orbit. While, since both u and v keep oscillating for the blue 
dotted curves in Figures 2(c)-(d), the blue dotted curve in Figure 2(b) is then a precessing 
elliptical orbit. 

Figure 3 is the case when we set M^, = 0.1 and k2 = 0.3. We assume the frictional 
parameter a = 0.01 for this case. From Figure 1, we know that there is only one fixed point 
and it is close to the inner edge of the belt, where u— 1/ri = 1/3. Figure 3(a) shows that 
this fixed point is a spiral point. The solution curve on the u — v phase-plane is then plotted 
as an orbit on the x — y plane in Figure 3(b). We can see that the planet moves across the 
belt initially and gradually settles down on a circular orbit close to the inner edge of the 
belt. Figures 3(c)-(d) show u and v as functions of 9 for this orbit. 

Figure 4 is the case when we set — 0.5 and k2 — 0.06. We assume the frictional 
parameter ck = 0.1 for this case. We know, from Figure 1, that there are two fixed points. 
Figure 4(a) shows that the one around m = 0.1 is a center point and another near the 
outer edge of the belt is a spiral point. The blue dotted curve around the center point in 
Figure 4(a) is then plotted on the x — y plane as the blue dotted orbit in Figure 4(b). It is 
obviously a precessing elliptical orbit. The red solid curve approaching the spiral point in 
Figure 4(a) is then plotted on the x — y plane as the red solid orbit in Figure 4(b). This 
orbit keeps changing the semimajor axis and eccentricity and finally becomes a circular 
orbit with radius r = r2, i.e., close to the outer edge of the belt. Figures 4(c)-(d) show u 
and V as functions of 9 for both above orbits. 

Figures 5-6 are for the case when we set = 0.1 and k2 — 0.1. We assume the 
frictional parameter a = 0.1 for this case. We know from Figure 1, that there is only one 
fixed point and Figure 5(a) shows that this fixed point is a center point. We found that 
there is a "hmit cycle" that all solution curves out of this cycle are approaching this cycle. 
However, all the solution curves within this cycle do not, they still behave as normal curves 
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around a center point. Therefore, this kind of "hmit cycle" is called a "semi-attractor" . 

To understand how these outer solution curves approach this "limit cycle" , we simply 
pick up the blue dotted curve's representative points which are on the right hand side of 
the fixed point and cross v = line, i.e. the points satisfying w > 0.12 and v = in Figure 
5(a). Figure 5(b) is the plot of the u values of these points, u, as a function of 9 when we 
set initial ^ = 0. We can see that u approaches tO'U = l/r2 = l/6 and the approaching 
speed decays very quickly. In Figure 5(b), we see that the curve almost becomes a straight 
line, i.e. du/d9 = when 9 > 100. 

The blue dotted curves in Figure 5(c)-(d) show u and v as functions of 9 and the 
curve in Figure 6(a) shows the orbit on the x — y plane for the above solution curve which 
approaches the "limit cycle". Figure 6(b) is the closer view of Figure 6(a). 

Moreover, the red solid curves in Figure 5(c)-(d) show the u and v as function of 9 and 
the curve in Figure 6(c) shows the orbit on the x — y plane for the red solid solution curve 
around the center point in Figure 5(a). 

Figure 6(d) is another orbit for a given solution curve approaching the "limit cycle". 
The orbit is getting closer to a precessing eUipse when the solution curve approaches the 
"limit cycle" as we can see in both Figure 6(b) and Figure 6(d). 

5. Conclusions 

We have studied the orbital evolution and also stability of planetary systems with 
interaction from belts by the standard phase-plane analysis. We regard the mass of the belt 
Mj, and k2 {— rri^/l^) as the two main parameters to determine the location of fixed points 
and also classify the outcome of orbital evolution. Please note that the fixed points on the 
phase space corresponds to the circular orbits on the orbital plane. 
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Prom the results in Figure 1, we know that when the belt's mass is small, the bifurcation 
diagram is almost a straight line, i.e. Keplerian orbits. When the belt's mass is larger, there 
could be more than one fixed point, one still corresponding to the Keplerian orbit and the 
new fixed points always close to the inner edge u — 1/ri or outer edge u — l/r2 oi the belt. 

However, since the needle-like structure around the outer edge in Figure 1 is very tiny, 
the probability that the planet moves stably around the outer edge is much smaller than 
near the inner edge. This conclusion is consistent with the principle result in Jiang & Yeh 
[2003]. 

There is one interesting case in which we found a limit cycle of semi-attractor type: the 
solution curves lying outside of the cycle approach this cycle but the ones within this cycle 
do not. In this case, the orbit is getting similar to a precessing ellipse when the solution 
curve approaches the limit cycle. 

What could we learn for the Solar System from these theoretical results ? 

It is known that there is a belt, the Kuipcr Belt, in the outer Solar System after the 
first Kuiper Belt Object (KBO) was detected (Jewitt & Luu [1993]). Many more KBOs 
have been detected since then and now there are more than 500 discovered KBOs. Allen, 
Bernstein & Malhotra [2001] did a survey and found that they could not find any KBOs 
larger than 160 km in diameter beyond 50 AU in the outer solar system. They suggested 
several possible explanations for this observational result. 

If we apply our model to this problem, the point mass which represents the planet in 
our model, now represents a KBO. Our results imply that the big KBOs would have greater 
probability of approaching the inner edge of the Kuiper Belt and the small KBOs would 
not. This is because the value of k2 would be larger when m is large for a given initial 
angular momentum, and thus it corresponds to the right branch of curves in Figure 1. If 
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m is much smaller and k2 < 0.1, it corresponds to the left branch of curves in Figure 1 and 
thus the small KBOs might not approach the inner edge of the Kuiper Belt. 

Therefore, our result is consistent with one of the possible reasons: larger KBOs might 
have been displaced to their present orbits by a large-scale rearrangement of orbit (Allen, 
Bernstein & Malhotra [2001]). In fact, our results provide a natural mechanism to do this 
orbit rearrangement: larger KBOs might have been moving towards the inner edge of the 
belt due to the friction from the belt. 
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A. Gravitational Force from the Belt 



Let F(^) be an elliptic integral of the first kind and assume 



(Al) 



By sin^ |(1 - cos 2A), we have 



Tr/2 



2 sin^ 6' 



I (r + r') n 



2 Jo _ 2ff' cos(0) 



#. (A2) 



Therefore, from Equation (A2), we have 

1 



27r 



_|_ J.I2 _ 2ff' cos(0) _|_ _ 27-7-' COS(0) 



(A3) 



The gravitational potential is 

r2Tv p(r')r' 



V{r) = -G / 

Jri Jo 



Iri Jo ^j.2 _^ J.I2 _ 2rr'cos(0') 
where Equation (A3) has been used. 



dr'dcf)' 



-AG / 

Jr\ 



F{i)p{ry 



r + r' 



dr', (A4) 



Now we differentiate this potential with respect to r to get the gravitational force, so 



h--E--^G ' dr' + AG FiOp{ry-^(^]dr\ (A5) 

or Jri r + r Jn or\r + rJ 



Let E{^) be an elliptic integral of the second kind and ^' = ^/l — Since F(^) is an 
elliptic integral of the first kind, from Byrd & Friedman [1971], we have 



dFiO _E-eF{0 _ E 



e(i - e) e 



(A6) 



By Equation (Al) and (A6), we calculate 



d^ dr 



E 



l^-f 7.' _ 7. E(r + r') Fij' — r') 



r (r + r')^ 2r(r — r') 2r(r + r') 



(A7) 
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We substitute Equation (A7) into Equation (A5), so we have 



= -— = -2G r 

dr Jri r 



dr'. 



(A8) 



If the planet is on the belt region, there is a singularity for and we introduce a small 
e > to avoid it. That is, we numerically calculate the force from the belt by the following 
formula: 



rr 



p{r')r' 



dr 



rr2 

' -2G 

Jr+( 



p{r')r' 



r+e r 



dr' . 
(A9) 



Similar treatment can be done for the integrals in Appendix B. 



B. Derivatives of 



From Equation (A8), we calculate 



dh{r) 
dr 



2G , 

/ Jri 

2G 

r Jr\ 



' E+^F 



dE dj , 

dr _^ 



dr' 
1 



dF rfg , 



dr'(Bl) 



From Byrd & Friedman (1971), 



dE 1,^ r + r' ^ 

dt. t. 2Vrr' 



Thus, 



dE d^ f' — f 



d^ dr 2r(r + r') 

We substitute Equation (A7) and Equation (B2) into Equation (Bl), so we have 



(B2) 
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Jri (r — r'j^(r + r'j r + r J 



9r 

Hence, from Equation (A8), we have 



(B3) 



ar u^ u W^^'' '{{l-r'uf{l + r'u) 1 + r'u j ^ ' 
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Fig. 1. — The bifurcation diagram of fixed points on the k2 — u plane with ri = 3, r2 
where the green dashed curves are for Mt = 0.01, the blue dotted curves are for Mb = 
the purple dotted curves are for Mb = 0.5 and the red solid curves are for Mb = 1. 
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Fig. 2. — The result for Mf, = 0.5 and k2 = 0.6 with a = 0.001, ri = 3, r2 = 6. (a) Solution 
curves on the u — v phase-plane, (b) Orbits on the x — y plane, (c) The plots of m as a 
function of 9. (d) The plots of f as a function of 9. The red solid (blue dotted) curves in all 
panels correspond to the same solution. 
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Fig. 3. — The result for Mb = 0.1 and k2 = 0.3 with a = 0.01, ri = 3, r2 = 6. (a) A solution 
curve on the u — v phase-plane, (b) An orbit on the x — y plane, (c) The plot of m as a 
function of 6. (d) The plot of f as a function of 6. The curves in all panels correspond to 
the same solution. 
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Fig. 4. — The result for Mf, = 0.5 and ^2 = 0.06 with a = 0.1, ri = 3, r2 = 6. (a) Solution 
curves on the u — v phase-plane, (b) Orbits on x — y plane, (c) The plots of m as a function 
of 6. (d) The plots of f as a function of 6. The red solid (blue dotted) curves in all panels 
correspond to the same solution. 
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Fig. 5. — The result for = 0.1 and k2 = 0.1 with a = 0.1, ri = 3, r2 = 6. (a) Solution 
curves on the u — v phase-plane, (b) -u as a function of 6, where u are the values of u of the 
blue dotted curve in (a) and satisfied m > 0.12 and v = 0. (c) The plots of u as a function 
of 6. (d) The plots of f as a function of 6. The red solid (blue dotted) curves in all panels 
correspond to the same solution. 
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Fig. 6. — The result for Mb = 0.1 and k2 = 0.1 with a = 0.1, ri = 3, r2 = 6. (a) The orbit 
on the X — y plane for the blue dotted curve in Figure 5(a). (b) The closer view of (a), (c) 
The orbit on the x — y plane for the red solid curve in Figure 5(a). (d) Another orbit for 
arbitrary given solution curve approaching the limit cycle. 



